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A  numerical  procedure  based  on  the  time-dependent  Schrddinger  equation  for  the  modeling 
of  multidimensional  mesoscopic  devices  is  presented.  The  primary  features  of  this 
numerical  method  are  a  tight-binding  formulation  of  the  quantum  mechanical  Hamiltonian, 
an  alternating  direction  implicit  generalization  of  the  Crank-Nicholson  method  for 
solving  the  discretized  multidimensional  time-dependent  Schrodinger  equation,  and  a 
numerical  implementation  of  absorbing  boundary  conditions  for  propagating  wavefunctions  at 
the  open  boundaries  of  the  simulated  region.  To  show  the  capabilities  of  the  absorbing 
boundary  scheme,  numerical  results  are  presented  for  the  diffusion  of  localized  wavcpackets 
out  of  open  simulation  regions.  As  an  application  to  mesoscopic  systems,  results  arc 
presented  for  the  switching  in  a  /'-structure  quantum  modulated  transistor  with  a  continuous 
input  vvavefunction. 


Aconfljiua  For 


MT'7  WAkl 
anc  TAH 
UmuunctfKojd 
Juotificat  lott. 


4 


Bv ... 

Distribution/ 

Availability  Cod' 

Diat 

(Avail  and /•* 

1  Spnol&l 

■ 

PvA  ! 

_ 

I.  INTRODUCTION 

The  transient  behavior  of  one-dimensional  quantum 
devices,  such  as  the  resonant  tunneling  diode,  often  is  sim¬ 
ulated  using  numerical  Wigner  function  methods.1-2  This 
technique  is  appealing  because  the  Wigner  distribution 
function  can  be  regarded  to  some  extent  as  the  quantum 
equivalent  of  a  classical  distribution  function.  However, 
since  both  spatial  coordinates  and  momentum  are  simulta¬ 
neously  considered,  for  the  simulation  of  an  //-dimensional 
structure  the  Wigner  function  is  2/i-dimensional.  Thus,  the 
simulation  of  novel  quantum  devices,  such  as  the  '/’-struc¬ 
ture'1  and  Aharonov-Dohm  effect4  transistors,  quickly  be¬ 
comes  prohibitive  using  the  Wigner  function  approach  be¬ 
cause  at  least  two  spatial  dimensions  must  be  included  in 
the  simplest  models.  An  alternative  approach  is  the  direct 
solution  of  the  //-dimensional  Schrddinger  equation.  While 
a  summation  over  incident  states  is  required,  this  summa¬ 
tion  will  be  limited  in  many  cases.  In  this  paper  a  numer¬ 
ical  method  suitable  for  l he  solution  of  the  multidimen¬ 
sional  time-dependent  Schrddinger  equation  in  mesoscopic 
structures  is  presented.  To  provide  an  unambiguous  treat¬ 
ment  of  quantum  mechanics  for  the  spatially  discretized 
coordinate  system  that  must  be  employed,  a  tight-binding 
formalism  is  used.  To  ensure  stability  and  unitary  time 
evolution  of  the  wavefunclion,  a  half-implicit  Crank- 
Nicholson  scheme  is  employed.  To  allow  simulation  of 
transport  through  device  structures,  absorbing  boundary 
conditions  are  included  within  this  scheme.  Results  are 
presented  first  for  the  diffusion  of  localized  wave-packets 
out  of  the  open  simulation  regions,  to  demonstrate  bound¬ 
ary  absorption  of  transient  current  flow.  Then  results  u re¬ 
presented  for  switching  in  a  /'-structure  quantum  modu¬ 
lated  transistor,  where  the  quantum  interference  puttern  of 
the  wave-function  is  varied  by  changing  the  effective  length 
of  a  side  arm  in  a  quantum  wire.  In  this  last  example,  in 
addition  to  boundary  absorption,  a  continuous  incoming 
vvavefunction  is  included  in  the  boundary  conditions. 


II.  NUMERICAL  APPROACH 
A.  Discretization  of  Schrodinger’s  equation 

The  time  evolution  of  the  spatially  and  temporally  con¬ 
tinuous  carrier  wavefunclion  i/»(r,t) s (r| i/»(t)>  is  emu¬ 
lated  by  that  of  a  discretized  wavefunclion  (R|  0(t„))  on  a 
uniform  virtual  tight-binding  lattice.  The  tight-binding 
Hamiltonian  used  here  is  of  the  form* 


//  X  I  R>c-k(R 

K 


(I) 


+•  X  <|R>v<R+JA|  f  | R)v(R— jA| ), 

where  R  are  the  lattice  sites,  A  is  the  lattice  spacing,  j  are 
the  set  of  orthogonal  unit  vectors  of  the  multidimensional 
space,  and  R±jA  are  the  corresponding  nearest  neighbor 
sites,  as  depicted  in  Fig.  1.  The  “hopping  potential"  v  and 
the  site  potentials  eH  arc  defined  by 

v.r  fi2/2/n*A2 


and 


(2) 


(3) 


eH  X  2  v+V{R), 

I 

respectively,  where  K(R)  is  the  potential  energy.  Thus,  by 
applying  this  tight-binding  Hamiltonian  to  the  wavefunc- 
tion  < R 1 1//) ,  with  the  condition  (R|R’)~=5Hih',  one  obtains 

<R[//f d>>  -  -~T>  X  (<R  +  JA|tf/>  -  2 < R 1 1//> 
2w*A‘  j 

I  <R  — jA|t//»  4-  F(R)<R|0).  (4) 

This  result  is  equivalent  to  a  uniform  mesh  finite  difference 
discretization  of  the  continuum  effective  mass  Humiltonian 
// -  {fi2/2m*)V2  +  F(r).  The  tight-binding  formalism  Is 
used  here  because  it  provides  an  unambiguous  treatment  of 
quuntum  mechanics  for  the  spatially  discretized  system 
that  must  bo  employed,  and  because  it  provides  computi¬ 
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F'lG.  I.  Simulation  region  of  the  virtual  tight-binding  lattice  showing  the 
coupling  of  sites  to  their  neighbors.  Departures  from  the  nearest-neighbor 
coupling,  Eq.  ( I ),  occur  at  the  boundaries,  particularly  for  the  absorbing 
boundaries  where  complicated  coupling  patterns  may  result  from  the 
extrapolation  schemes  of  Eqs.  ( !l ) — f  10). 


bility  with  a  previously  developed  Green’s  function  method 
for  calculating  steady-state  transmission  coefficients.6  As 
will  he  shown  later,  inconsistent  approximations  of  the 
continuum  are  a  potential  source  of  error  in  absorbing 
boundary  conditions. 

It  is  known  that  a  Crank-Nicholson  approach  should 
he  used  for  stability  and  preservation  of  the  wavefunction 
norm,  when  a  first  order  forward  difference  is  used  for  the 
time  discretization  of  the  Schriklingcr’s  equation.7  The 
Hamiltonian  is  approximated  by  the  average  between  the 
eurrcait  and  the  next  time-step  value,  yielding  the  half- 
implicit  discretized  equation 

[!/'(•/„+  At) )-  |  0(/„ ) )  U\W„) >  +  H \Mt„  +  At)) 

1  kt  "  2 

(5) 


where  A t  is  the  time  step.  For  multidimensional  calcula¬ 
tion,  however,  large  matrices  are  generated.  To  avoid  this 
complication,  the  space  operator  can  be  split  into  compo¬ 
nents  and  an  alternating  direction  implicit  (ADI)  ap- 
proachl',<)  employed.  Specifically,  in  two-dimensions  (2-D), 
two  subiterations  with  time-step  At/2  are  used, 

■)-  vAr )  >  +  At ) ) 

+  //,!</'(/,))] 

if  |  1 1>(  t„  +  A/  )  >  r. ,  Hi  i  M  in  -I-  vA  /  )  )  -1-  JA/  t  II,  I M  tn  +  'At)) 
+  lly\  -I  A/)>],  (6) 

where  flx  and  Hy  result  from  the  decomposition  of  the 
tight-binding  Hamiltonian  of  Eq.  (1)  into  its  orthogonal 
components.  The  structure  of  the  scheme  is  given  in  (6) 


*'te  .v-  anc’ f-cr.mpotieitts  are  a'lernaf ive'y  im- 
>i  •  I  in  ihe  two  subiterations.  This  procedure  require*  only 
I  .  '  ut ion  of  tridiagonal  matrices  along  the  horizontal  or 
vertical  niedt  lines.  Furthermore,  since  such  sets  of  matrix 
■a  mi  nis  arc  decoupled  from  each  other,  a  very  efficient 
itnplerne  :t;:t  oti  on  a  supercomputer  is  possible. 


13  Absorbing  boundary  conditions 

Ileal  i*rc  modeling  of  mesoscopic  devices  requires  sim¬ 
ulation  of  both  open  and  closed  boundaries.  Numerical 
simulation  o  ' closed  (reflecting)  boundaries  is  straightfor¬ 
ward;  the  eojpling  potential  v  to  the  boundary  tight-bind¬ 
ing  lattice  files  simply  can  be  reset  to  zero,  or  equivalently 
ihe  boundary  wavefunctions  fixed  at  zero.  The  challenge 
here.  ,i  challenge  common  to  all  fields  concerned  with  wave 
propagation, 10  is  numerical  simulation  of  outflow  of  wave- 
functions  through  the  simulation  region  open  boundaries 
without  introducing  spurious  reflections  due  to  the  numer¬ 
ical  scheme.  The  work  of  Mains  and  Haddad  in  one-di¬ 
mension  (1-D)11,12  suggests  that  a  curve-fitting  approach 
for  the  wavefunction  near  the  open  boundaries  may  pro¬ 
vide  a  solution  to  this  problem.  They  apptoximate  the 
boundary  wavefunctions  as 

t i> ,  n  (•*)=;  aeikx  -t-  (  b,+  c/x )  c  ~  kx 


tl>,su(x)--U>r  +  c,x)ek'*  (7) 

and  ther  apply  the  continuum  Hamiltonian  to  calculate 
the  time  evolution  of  the  wavefunction  at  the  boundary. 
Here,  k  is  the  wavevector  of  a  left  incoming  monochro¬ 
matic  wave,  k'  is  the  corresponding  wavevector  at  the. 
right-hand  side  boundary  with  an  allowance  for  the  poten¬ 
tial  energy  difference  between  the  boundaries,  and  the  co¬ 
efficients  b  and  c  are  obtained  from  a  linear  fit  to  the  wave- 
function  at  the  boundary  and  first  internal  points  of  the 
discrete  simulation  grid  employed.  This  method  was  effec¬ 
tive  for  their  1-D  simulations  of  resonant  tunneling  struc¬ 
tures.  However,  this  boundary  condition  loses  accuracy  as 
the  relative  variation  in  the  wavefunction  momentum  from 
the  specified  monochromatic  limit  increases.1'1  Further¬ 
more,  the  application  of  the  continuum  Hamiltonian  to  the 
boundary  wavefunctions  of  Eq.  (7)  is  inconsistent  with  the 
finite  difference  discretization  of  the  Hamiltonian  used 
within  the  domain.  A  discrepancy  exists  between  the  par¬ 
abolic  band  structure  of  the  continuum  Hamiltonian  and 
the  tight-binding  band  structure  implicit  to  the  finite  dif¬ 
ference  discretization.  This  discrepancy  is  particularly  sig¬ 
nificant  for  large  wavevectors. 

In  this  paper  a  more  general  and  robust  approach  to 
fitting  the  boundary  wavefunction  and  evaluating  the 
boundary  Hamiltonian  is  presented.  The  wavefunction  is 
extrapolated  on  the  tight-binding  lattice  sites  immediately 
outside  the  open  boundaries  using 


<Rr I  </>(!))  =*  T[<R|  0( /„))]}' 

or,  similarly, 


<R|*(/» 

<R|*(f„)> 


<Re[ t/'(0> » P[(R| 0(t„)> ]  -f  y((R|0(f)> 


(8) 
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(9)  probability  How  from  an  external  lattice  site  to  a  boundary 
site  has  been  restricted.  It  can  be  shown  readily  that 


where  Rc  is  an  external  lattice  site,  and  F  and  y  are  inde¬ 
pendent  extrapolation  functions-.  For  the  2-D  ADI  scheme 
(b),  i  r,t  |  ‘A/  or  /  —  i„  -f-  A /  =  t„  ,  The  tight-binding 
Mamijloninn  of  Eqs.  ( I  )-(4)  then  is  applied  to  the  bound¬ 
ary  sites,  just  as  it  is  applied  to  the  internal  sites,  to  obtain 
the  time  evolution  of  the  boundary  wavefunction,  as  de¬ 
picted  in  Fig.  I.  Inclusion  of  specified  incoming  wavefune- 
lions  (R|t/’;}  in  the  boundary  conditions  is  accomplished 
by  replacing  < R((.,  | «//)  in  Eqs.  (8)  and  (9)  with 
<R(,  l|i/’)  < R(,.)  i  ifri).  Further,  both  of  the  extrapolation 

schemes  I’  and  y  can  be  varied  as  a  function  of  R,,  and  /„  if 
required. 

T  he  separation  of  the  extrapolation  schemes  of  Eqs. 
(8)  and  (9)  into  the  components  P  (explicit)  and  y  (im¬ 
plicit)  is  intended  to  satisfy  two  conflicting  requirements: 
Accuracy  of  the  extrapolation  and  realizability  of  the  im¬ 
plicit  scheme  (6).  Because  the  primary  extrapolation  T  is 
performed  Only  explicitly,  complicated  extrapolation 
schemes  can  by*  used  without  reducing  the  efficiency  of  the 
ADI  numerical  approach  described  in  Sec.  II  A  above. 
This  flexibility  allows  complicated  admixtures  of  propagat¬ 
ing  and  evanescent  states  to  be  fit  at  the  boundaries.1415  On 
the  other  hand,  )'  is  required  to  extrapolate  only  small 
changes  in  the  boundary  wavefunction  over  the  time  step 
and  thus,  extrapolation  schemes  that  can  be  readily  evalu¬ 
ated  implicitly  are  sufficient.  Also,  here  there  is  no  inherent 
restriction  on  the  nature  of  the  incoming  wavefunction; 
<  R  1 1/',)  can  be  monochromatic  or  otherwise,  c.g.,  a  prop¬ 
agating  wavepaeket.  In  contrast,  in  Refs.  II  and  12  the 
incoming  wavefunction  is  necessarily  monochromatic.  Fi¬ 
nally.  application  of  the  tight-binding  Hamiltonian  of  Eqs. 
( 1 )— f  4 )  to  the  boundary  lattice  sites  as  to  the  internal  sites 
allows  no  possibility  of  discrepancy  in  the  resulting  energy 
band  structure  between  the  two  sets  of  lattice  sites. 

For  the  example  calculations  presented  in  the  next  sec¬ 
tion,  \ve  chose  the  following  simple  scheme,  in  the  form  of 
Eq.  (8): 


( R /,  I  /<  A  |  ~ 


n  <r„  /!  «a  i  </'(/„)>"« 

n  n 

V  <R/,  --  fin  A 1 1/;(/)) 

L  <>„ 

n  0 


<R/,-  /t/i A  |  (/’(/„)) 


(10) 


where  R,,  is  the  boundary  site  and  fi  is  a  unit  vector  normal 
to  the  boundary  and  directed  outwards,  as  shown  in  Fig.  I. 
An  appropriate  choice  of  the  coefficients  a„  provides  any 
desired  polynomial  curve  fit  near  the  boundary.  For  exam¬ 
ple,  for  a  linear  curve  fit,  of)  —  2,  0\  -  -  I,  an  —  0,  and 
for  a  parabolic  curve  fit,  un  —  i,  c/|  —  —3,  —  I, 

</„  ,  0.  In  addition,  to  prevent  the  "absorbing”  bound¬ 

aries  from  inadvertently  becoming  probability  sources. 


^ |  <  R |  ■/')  | 2  -~jr  X  lm[<0|R><R  t  JA|  V'> 

f  <«/' |R><R . j A | V'> ]  (ID 

for  the  tight-binding  formalism  used  here.  Therefore, 
should  the  phase  of  P((R| 4>(t„))  —  (R|  V'/(f„)))  ever  lag 
that  of  <R*f  «//(/„))  —  (R/,|0/(7„))  for  any  boundary  lattice 
site,  we  project  the  former  complex  quantity  onto  fhe  latter 
to  obtain  a  new  prediction  T  of  equal  phase,  with  the  min¬ 
imum  change  to  the  initial  prediction  possible. 


III.  NUMERICAL  EXAMPLES 

A.  Diffusion  through  absorbing  boundaries 

The  quantum  mechanical  diffusion  of  localized  wave- 
packets  (nit  of  open  regions  has  been  simulated  to  demon¬ 
strate  the  effectiveness  of  the  absorbing  boundary  scheme 
presented  here.  Because  of  the  wide  spectral  spread  in  en¬ 
ergies,  small  mean  momentums,  and  open  boundaries  on 
all  sides,  these  diffusion  examples  pose  in  many  respects  a 
difficult  numerical  challenge.  For  the  first  such  example, 
the  diffusions  of  two  identical  stationary  Gaussian  wave- 
packets  out  of  concentric  square  domains  were  simulated 
and  the  results  compared,  as  shown  in  Figs.  2(a)  and  2(b). 
Here,  the  probability  distribution  j  <R  1 1/')  | 2  obtained  for 
the  smaller  domain  is  overlaid  on  that  obtained  for  the 
larger  to  facilitate  comparison.  The  difference  between  the 
two  results  at  the  boundary  of  the  inner  region  is  a  measure 
of  the  error  induced  by  the  application  of  absorbing  bound¬ 
ary  conditions  of  Eq.  ( 10)  with,  in  this  case,  a  cubic  curve 
fit.  After  200  fs  of  simulation  such  error  remained  negligi¬ 
ble,  as  to  be  expected;  for  the  primary  extrapolation  T,  a 
parabolic  or  higher  order  curve  fit  is  exact  for  a  Gaussian 
wavepaeket.  In  this  and  all  following  examples  the  effective 
mass  in*  of  GaAs  (0.067m,,),  a  light-binding  grid  spacing 
A  of  If)  A  and  a  time-step  At  of  1.0  fs  were  used.  For 
comparison,  for  the  wavefunction  of  this  example,  which 
has  zero  mean  momentum  (and  no  reference  incoming 
wavefunction),  the  absorbing  boundary  condition  in  Refs. 
11  and  12  becomes,  in  effect,  a  reflecting  condition  as  the 
continuum  Hamil'onian  vanishes  at  the  boundary,  allow¬ 
ing  no  time  evolution  of  the  wavefunction  there. 

For  the  second  example,  the  diffusions  of  identical 
wavefuuctions,  each  this  time  initially  a  sum  of  two  Gauss¬ 
ian  wavepackets  with  a  relative  phase  difference  of  it/2 
radians,  out  of  two  concentric  rectangular  domains  were 
simulated  and  compared,  as  shown  in  Figs,  2  (e)  and  2(d). 
In  these  simulations,  the  curve  fit  of  Eq.  (10)  wus  no 
longer  exact,  and  one  can  notice  small  discrepancies  be¬ 
tween  the  two  solutions  at  the  inner  boundary  after  200  fs 
of  simulation.  The  error  exhibited  in  Fig.  2(d),  however,  is 
emphasized  by  the  choice  of  scale  which  is  50  times  smaller 
than  in  Fig.  2(c). 


7155 


J.  Appl.  Phys.,  Vol  69,  No.  10,  15  May  1991 


Register,  Ravatolt,  and  Hess  7155 


I'lO,  2.  Simulation  n-sulls  lor  llie  diffusions  of  Gaussian  wuvcpackels  out 
of  concentric  open  regions:  (a)  initial  probability  distribution  and  (b) 
probability  distribution  after  2(X)  fs,  for  a  single  stationary  Oaussiun 
wave-packet;  and  f c>  initial  probability  distribution,  and  (d)  probability 
distribution  after  200  fs,  for  a  wavefunction  obtained  from  the  sum  of  two 
Gaussian  wavcpackcls  with  a  relative  phase  difference  of  it/ 2  rudians. 
The  boundaries  of  the  inner  regions  are  outlined  when  not  otherwise 
apparent 
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B.  Switching  of  a  quantum  Interference  structure 

As  a  final  example,  the  2-D  numerical  scheme  pre¬ 
sented  here  was  used  to  simulate  the  switching  transient  in 
a  7-structure  quantum  modulated  transistor.3  An  idealized 
2-D  structure  delimited  by  hard  walls  was  assumed  to  al¬ 
low  comparison  of  the  results  obtained  here  with  steady- 
state  transmission  coefficients  obtained  using  the  steady- 
state  Greens  function  method  of  Ref.  6.  However,  this 
assumption  does  not  reflect  a  limitation  of  the  method; 
finite  potential  walls  could  have  been  modeled  via  the  spa¬ 
tially  varying  site  potentials  of  Eq.  (3).  The  input  and 
output  leads  of  the  structure  were  modeled  as  quantum 
wires  with  a  width  W  of  100  A.  The  side  arm  had  the  same 
width  and  an  initial  length  L  of  70  A.  Again,  the  boundary 
conditions  of  Eq.  (10)  were  used,  though  this  time  a  sim¬ 
ple  linear  curve  fit  proved  sufficient.  However,  an  incoming 
wavefunction  (R  |  t/r,)  <x%(  Y)exp(ikX)  was  included  in  the 
left  boundary  condition.  Here,  Y)  >s  the  ground  state 
transverse  wavefunction  of  the  lead  and  k  is  the  incident 
momentum  wavevector  taken  to  be  27T/200  A  ~  '.  The  ini¬ 
tial  steady-state  probability  distribution  |<R|^)|2  is 
shown  in  Fig.  3(a).  Due  to  quantum  interference  effects  in 
the  stub  region,  the  transmission  probability  through  the 
structure  was  nearly  zero.  One  can  see  a  pronounced  peak 
of  probability  inside  the  stub  region,  and  a  standing  wave 
pattern  in  the  input  lead  indicating  nearly  total  reflection. 
The  time  evolution  of  the  probability  distribution  after 
switching  the  stub  length  from  70  to  100  A  is  shown  in 
Figs.  3(b)— 3(i).  Here,  a  uniform  probability  scale  is  main¬ 
tained  throughout.  The  transient  exhibits  a  quick  collapse 
of  the  probability  peak  in  the  stub  and  a  corresponding 
wavefront  in  the  output  lead.  The  decrease  of  the  reflected 
wave  in  the  input  lead  is  evident  from  the  progressive  re¬ 
duction  of  the  standing  wave  pattern.  By  500  fs  the  struc¬ 
ture  was  clearly  on  with  a  transmission  probability  near 
unity,  with  only  a  hint  of  reflection  at  the  input  noticeable 
from  the  weak  standing  wave  pattern.  Only  a  small  frac¬ 
tion  of  the  original  probability  peak  inside  the  stub  re¬ 
mained.  Beyond  this  time,  although  switching  was  com¬ 
plete,  the  time  evolution  still  exhibited  small  persistent 
oscillations.  The  abrupt  initial  switching  excited  localized 
eigenstates  of  the  stub  at  various  energies.  These  eigen¬ 
states,  which  are  evanescent  in  the  leads,  did  not  decay  in 
time  since  dissipation  is  not  included  in  the  model,  and 
thus  continued  to  beat  against  the  propagating  state.  If 
desired,  such  localized  states  can  be  filtered  out  by  Fourier 
transforming  from  time  to  energy  at  the  energy  of  the  in¬ 
cident  wave.  Indeed,  this  approach  was  used  to  obtain  the 
initials  conditions  for  this  simulation,  Fig.  3(a),  after 
switching  from  a  simple  quantum  wire  with  no  stub.  On 
the  other  hand,  localized  states  can  be  isolated  by  Fourier 
transforming  at  their  energies.  The  simulation  results  were 
not  appreciably  affected  by  changing  the  lengths  of  the 
leads  to  check  the  absorbing  boundary  conditions.  A  sim¬ 
ilar  simulation  of  switch-off  exhibited  a  somewhat  quicker 
and  less  oscillatory  switching  characteristic.  In  both  cases, 
the  steady-state  transmission  probabilities  obtained  here 
are  identical  to  those  obtained  using  the  steady-state 
Green's  function  method  of  Ref.  6. 
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IV.  CONCLUSIONS  absorbing  boundary  conditions,  suitable  for  the  simulation 

We  have  presented  a  numerical  method  based  on  the  of  transient  response  in  multidimensional  mesoscopic  dc- 
time-dependent  Schrddingcr  equation  and  incorporating  vices.  We  have  shown  results  for  the  diffusion  of  localized 
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wavcpae'iets  out  of  open  syslems  and  for  the  switching  of 
a  /-structure  quantum  modulated  fransistor.  Statistical 
distributions  in  the  state  occupancy  in  the  contact  reser¬ 
voirs  can  be  included  by  a  weighted  sum  of  results  for 
monochromatic  incident  states  over  the  appropriate  energy 
range.  Future  work  includes  the  introduction  of  potential 
energy  self-consistency,  via  the  solution  of  Poisson’s  equa¬ 
tion,  and  of  a  spatially  varying  effective  mass. 
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